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Final  Report 


1.0  Introduction 

The  objective  of  Phase  1  of  this  project  was  to  demonstrate  feasibility  of  exploiting  the 
commonality  in  the  structure  of  the  governing  equations  in  fluid  dynamics  and  electromagnetics 
to  develop  a  higher-order  Magneto  Gas  Dynamic  (MGD)  solver  based  on  the  3D,  discontinuous 
Galerkin  (DG),  CEM  (Computational  ElectroMagnetics)  code  TEMPUS  [1]  under  development 
at  HyPerComp.  The  MGD  solver  was  to  be  restricted  to  inviscid  flows.  With  a  view  to 
achieving  this  objective,  we  first  investigated  implementation  of  a  DG  scheme  for  solving  the 
one-dimensional  (ID)  Euler  equations,  the  governing  equations  for  the  flow  of  an  inviscid  fluid 
in  ID.  The  Euler  equations  represent  one  of  the  simplest  non-linear  set  of  partial  differential 
equations  (PDE)  of  relevance  to  fluid  mechanics  and  wave  propagation.  We  studied  application 
of  the  DG  algorithm  for  numerical  solution  of  ID  Euler  equations  as  applied  to  simulation  of 
flow  in  a  shock  tube,  flow  through  a  quasi- ID  convergent-divergent  nozzle,  and  propagation  of  a 
periodic  density  wave.  The  implementation  employed  was  based  on  a  novel  linearization 
technique  for  the  flux  terms  that  simplified  the  development  of  explicit  schemes  as  well  as  a 
point-implicit  scheme.  Based  on  the  linearization,  an  implicit  formulation  for  the  ID  Euler 
equations  was  developed.  The  DG  formulation  tested  both  monomials  and  orthogonal 
polynomials  as  basis  functions.  The  formulation  was  extended  for  inviscid,  ID,  MGD  equations 
and  solutions  were  obtained  for  the  Brio-Wu  shock-tube  problem  [2]  and  the  Wu  intermediate 
shock  problem  [3].  This  was  followed  by  the  conversion  of  the  3D  DG  CEM  code,  TEMPUS,  to 
a  perfect  gas  Euler  solver.  The  experience  gained  from  the  ID  studies  enabled  efficient  code 
conversion.  Inviscid  flow  over  a  2D  cylinder  was  used  as  the  test  problem  to  verify  the  veracity 
of  the  conversion.  Future  work  will  add  the  additional  equations  required  to  complete  the  MGD 
equation  set.  The  code  conversion  task  validated  our  initial  assumption  that  the  commonality  in 
the  structure  of  the  governing  equations  cast  in  a  conservation  form  in  many  disciplines  such  as 
CEM  and  CFD  may  be  exploited  in  developing  a  common  platform  for  solving  problems  in 
multi-disciplines.  Design  strategies  for  the  development  of  a  common  platform  for  multi-physics 
computation  were  also  investigated. 

In  the  following  sections  we  present  the  details  of  our  work  in  Phase  I .  The  fundamentals  of  a 
DG  scheme  are  described  in  section  2.  Development  of  a  DG  scheme  using  orthogonal 
polynomials  for  a  system  of  conservation  laws  in  ID  is  detailed  in  section  3.  Numerical  schemes 
developed  using  the  DG  scheme  for  simulation  of  inviscid,  ID  and  quasi- ID  flows  are  presented 
in  section  4.  Results  obtained  using  the  DG  scheme  for  some  classical  problems  are  discussed  in 
section  5.  An  implicit  formulation  of  the  DG  scheme  is  presented  in  section  6.  The  3D  DG 
Euler  equations  are  described  in  section  7.  Evaluation  of  some  integrals  involved  in  the  3D 
formulation  is  discussed  in  section  8.  Steps  involved  in  the  conversion  of  the  CEM  code  to  an 
Euler  solver  are  detailed  in  section  9.  Some  3D  test  results  are  shown  in  section  10.  Design  of  a 
common  platform  for  multi-physics  computation  is  discussed  in  section  1 1 .  Suggestions  for  the 
extension  of  the  work  done  under  Phase  I  are  presented  in  section  12. 
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2.0  Discontinuous  Galerkin  Algorithm 


In  this  section,  a  brief  description  of  the  discontinuous  Galerkin  algorithm  is  presented  by 
applying  it  to  the  following  generic  system  of  conservation  laws: 

^  +  V-F{Q)  =  S{x,y,z,t)  (2.1) 

ot 

where  Q  is  the  vector  of  conserved  quantities,  F  is  the  flux  tensor,  and  S  is  the  vector  of  source 
terms. 

The  computational  domain  (D)  is  partitioned  into  Ne  non-overlapping  subdomains  (elements) 
De.  Within  each  element,  the  solution  Q  is  projected  onto  a  local  basis  vector  Vg(x,y,z); 
e=0, . . . ,Ne- 1  (e  —>  element),  and  i=0, 1 , . . . ,n- 1 .  It  is  represented  by 

n-l  ^  ■ 

Q(x,y,z,t)=  Y,  Ql(t)Ve(x,y,z),  for {x,y,z)z  De  (2.2) 

/=0 


where  Q^it)  are  time-dependent  coefficients.  The  ‘Discontinuous’  in  ‘Discontinuous  Galerkin’ 

refers  to  the  property  of  eqn.(2.2)  by  which  the  approximation  for  conserved  quantities,  while 
continuous  within  an  element,  is  not  necessarily  continuous  across  neighboring  elements. 

The  unknown  coefficients  Qg{t)  are  computed  via  a  weak  formulation  of  eqn  (2.1),  obtained  by 

multiplying  it  with  each  of  the  basis  functions  Vg  (x,  y,  z)  and  then  integrating  by  parts  over  each 
element  De  to  obtain 
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De 


i  v7  i  r 

Ot 


w+  1 

8De 


WeF 


dA=  jSvgdV;  e=0,...,Ne-l,  i=l,2,...,n-l 


(2.3) 


De 


where  dDe  is  the  boundary  of  element  De.  The  ‘Galerkin’  in  ‘Discontinuous  Galerkin’  refers  to 
satisfying  a  weak  form  of  eqn.(2.1)  rather  than  eqn.(2.1)  itself. 

The  volume  and  boundary  integrals  in  eqn.(2.3)  may  be  evaluated  using  either  numerical 
quadratures  or  exact  integrals  depending  on  the  basis  functions  employed.  In  either  case  the 
boundary  integral  involves  evaluation  of  the  flux  tensor  F  at  points  on  the  boundary.  Since  Q  is 
not  necessarily  continuous  across  an  element  boundary,  we  need  a  procedure  for  computing  a 
unique  value  for  Fat  a  boundary  point.  While  F  can  be  evaluated  in  a  number  of  ways, 
computing  it  by  solving  the  Riemann  problem  based  on  the  ‘left’  and  ‘right’  states  for  Q 
computed  at  an  element  interface  using  eqn.(2.2)  for  the  ‘left’  and  ‘right’  elements  respectively 
has  produced  the  highest  phase  accuracy  in  our  numerical  experiments.  Solving  the  Riemann 
problem  on  the  boundary  of  an  element  is  physically  equivalent  to  computing  F  by  considering 
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only  those  waves,  whieh  are  propagating  towards  the  boundary,  and  ignoring  those  whieh  are 
propagating  away  from  it. 

One  of  the  most  popular  Riemann  solvers  employed  in  CFD  (eomputational  fluid  dynamics)  is 
the  approximate  Riemann  solver  developed  by  Roe  [4],  Computation  of  F  at  a  boundary  point 
(BP)  using  Roe’s  Riemann  solver  starts  by  computing  Q  at  the  point  from  eqn.(2.2)  for  the  two 

parent  elements  of  the  boundary  point.  Let  these  Q  values  be  denoted  by  Qi  and  Qji  where 
the  subscripts  L  and  R  refer  to  the  “Left”  and  “Right”  element  respectively.  The  flux  tensor  at 
the  boundary  point  may  now  be  written  as 

F(BP)='^[f(ql)+f{qk\-^\J\(qk -Qi)  (2.4) 


J  is  given  by 


J 


=  R 


A 


(2.5) 


IS  a 


-  -  dF  - 

where  R  and  L  are  the  “right”  and  “left”  eigen  vectors  of  the  Jacobian  J  =  and  A 

dQ 

diagonal  matrix  whose  elements  are  the  absolute  values  of  the  eigen  values  of  J.  Equation  (2.4) 
may  be  interpreted  as  computing  the  flux  vector  at  a  boundary  point  using  the  central  difference 

operator,  ^[F(g2,)+^(0i?)]  and  the  dissipation  term  ^J^Qr-Ql)-  Computation  of  the 

dissipation  term  involves  evaluation  of  the  eigen  values  and  eigen  vectors  of  the  Jacobian  matrix. 
In  order  to  avoid  expensive  evaluation  of  eigen  vectors  without  severely  compromising  the 
accuracy  of  the  scheme,  we  employ  a  simpler  form  of  eqn.  (2.4)  wherein  the  dissipation  term  is 

computed  as^lA-j^axll^i?  ~Ql)-  This  leads  to  the  following  expression  for  the  flux  vector  at  a 

point  on  an  element  boundary: 


FiBF)  =  I [f(Ql  )+  )]- W ife  -  Ql ) 


(2.6) 


We  refer  to  such  a  scheme  as  the  “scalar  dissipation”  scheme  because  the  computation  of  the 
dissipation  term  involves  only  a  scalar  coefficient  and  not  a  vector  coefficient  as  in  the  case  of 
Roe’s  Riemann  solver. 


The  choice  of  the  basis  functions  Vg(x,y,z); i=0,l,...,n-l  determines  a  particular  type  and 

accuracy  of  the  discontinuous  Galerkin  method.  One  can  choose  monomials  such  as  1 ,  x,  y,  z, 
X  ,  y  ,  etc.  or  other  functions  such  as  orthogonal  polynomials  as  basis  functions.  Obviously,  with 
a  larger  basis  function  set,  one  can  approximate  the  solution  more  accurately  within  each 
element,  but  at  the  expense  of  higher  computational  cost.  We  refer  to  the  method  as  DG-pn, 
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where  n  is  the  order  of  the  polynomial  representing  the  solution.  It  has  been  shown  in  the 
literature  that  a  DG-pn  scheme  will  be  of  order  at  least  (n+1),  and  of  order  up  to  (n+2)  in  certain 
cases. 

A  second  order  TVD  scheme  also  approximates  the  variation  of  Q  within  an  element  by  a  linear 
function  and  thus  has  nominally  the  same  order  of  accuracy  as  the  DG-pl  scheme.  The  main 
difference  between  the  two  schemes  arises  from  the  manner  in  which  the  gradient  terms 

(Gx  ’  Qy  ’  Qz  )  are  computed.  While  a  TVD  scheme  approximates  the  gradient  terms  using 

numerical  differencing,  a  DG-pl  evaluates  them  by  solving  a  set  of  PDEs  derived  using 
eqn.(2.3). 

While  TVD  schemes  are  limited  to  2"^*  order  accuracy,  DG  schemes  of  arbitrary  order  of 
accuracy  have  been  formulated.  Third  and  higher-order  ENO  (Essentially  non-Oscillatory) 
schemes  which  satisfy  the  TVD  property  almost  everywhere  have  enjoyed  only  very  limited 
success.  On  the  other  hand,  in  CEM  solutions  have  been  obtained  for  some  three  dimensional 
(3D)  configurations  using  DG  schemes  of  order  up  to  eleven  and  it  has  been  shown  that  it  is 
possible  to  substantially  increase  accuracy  and  decrease  tum-around  time  using  higher  order  DG 
schemes. 

3.0  DG  Formulation  for  ID  Conservation  Laws  using  Orthogonal  Polynomials 

In  the  initial  stages  of  the  contract,  formulation  for  ID  conservation  laws  using  both  monomial 
and  orthogonal  basis  functions  were  studied.  The  orthogonal  basis  functions  are  preferred  over 
the  monomials  since  they  present  two  major  advantages  over  the  monomials,  namely 

1 .  The  equations  for  the  coefficients  are  uncoupled. 

2.  The  formulation  is  hierarchical  in  the  sense  that  to  increase  the  order  of  the  formulation  from 
n  to  n+1,  one  needs  to  only  add  the  equation  for  the  coefficient.  This  property  follows 
from  the  fact  that  the  coefficients  are  uncoupled.  On  the  other  hand,  in  the  case  of 
monomials  addition  of  a  basis  function  modifies  the  governing  equations  for  all  the 
coefficients  necessitating  the  reformulation  of  the  scheme. 

Another  difference  between  the  use  of  monomials  and  orthogonal  polynomials  arises  from  the 
fact  that  while  in  the  case  of  monomials  Qg  corresponds  to  Q  at  the  element  center,  in  the  case 
of  orthogonal  polynomials  it  corresponds  to  the  element  average.  In  other  words,  in  the  case  of 
monomials  Q  is  expanded  about  the  element  center  while  in  the  case  of  orthogonal  polynomials 
the  expansion  is  about  the  element  average. 

The  orthogonal  polynomial  basis  functions,  v*  ),  employed  in  the  formulation  are  the  Legendre 
polynomials.  They  are  obtained  sequentially  starting  with  the  zeroth  order  polynomial  by 
satisfying  the  following  conditions: 
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(3.1) 


1 

-1 


it 

2/  +  1 


v'(l)  =  l 

(3.2) 

v'(^)  =  (-l)'v'(-^) 

(3.3) 

Eqns.  (3.1),  (3.2),  and  (3.3)  lead  to 

2  3 

v^{^)  =  l,  =  =  v^fe)=^^'~^,....,andsoon  (3.4) 

Following  the  presentation  in  section  2,  we  first  expand  the  dependent  variable  vector,  Q ,  within 
an  element  as 

Q{x)=  zeivife)  (3.6) 

1=0 

where 

^=2-^^  (3.7) 

h 

where  Xg  is  the  center  of  the  element  and  h  is  its  length.  The  dependent  variables  for  the 

numerical  scheme  now  consist  of  the  {n+1)  vectors  Qg,i  =  0,n.  The  weak  form  of  the 
governing  equations,  eqn.(2.3),  may  now  be  written  as 

+  (3.8) 

2i  +  \  ot  dq 

Equation  (2.6)  is  employed  in  the  computation  of  F(i)  and  F{-  l)  that  involve  computation  of 
flux  vector  at  an  element  boundary. 

Orthogonality  of  the  basis  functions,  results  in  decoupling  of  the  gg' son  the  left-hand  side 

of  eqn.  (3.8).  Orthogonality  also  simplifies  higher-order  implementations  in  the  sense  that  the 

^  M  _1_  1 

order  of  accuracy  may  be  increased  from  n  to  n+1  by  just  adding  the  equation  forgg  .  The 
equations  for  Qi ,  i=0,n,  remain  unchanged.  The  integrals  on  the  right-hand  side  (RHS)  of  eqn. 
(3.5)  may  be  evaluated  either  using  a  Gaussian  quadrature  or  by  expanding  the  flux- vector  F  and 
the  source-vector  S  about  Qg  .  The  Taylor  series  expansion  for  F  about  Qg  may  be  written  as 
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0, ,  ^  dfi  ^  d'^fi  ^qj -^qk  ,  d^fi  ^qj -Mk -^qi 


y;-(0=y;-(0.)+z^A^y+ii 


dq 


+  ZZZ- 


J-'HJ  j  k  -yj 


kdqjdqk  2!  j  i  dqjdqkdqi 


3! 


+ 


(3.9) 


where  fi  and  qi  refer  to  the  component  of  the  flux  vector  and  dependent  variable  vector 
respectively,  and 


^qi  =qi-qei  =  'L  qii'H^) 

7=1 


(3.10) 


qg  ■  being  the  component  of  . 


Using  eqns.  (3.9)  and  (3.10)  the  flux  vector  Fmay  be  expanded  in  terms  of  the  basis  functions 
as 


^fe)=  ZFii^e.Ql . .elj-vife)  (3.11) 

Substitution  of  eqn.  (3.11)  in  eqn.  (3.8)  reduces  the  evaluation  of  the  flux-vector  integral  on  the 


RHS  of  eqn.  (3.8)  to  the  evaluation  of  the  integrals 

1  dv^ 

j  ^  i=0,  n,  j=0,  n  (3.12) 

The  source-vector  integral  may  also  be  obtained  using  the  same  procedure. 

The  final  form  of  the  DG  equations  up  to  p3  are  the  following: 

i=0:  /i0O, +F(1)-F(-1)  =  2-S0  (3.13) 

i=l:  !Iq\^+F{\)  +  F{-\)  =  2-F0  +  ^-S\  (3.14) 

i=2:  ^02, +F(l)-F(-l)  =  2-Fl  +  ^-52  (3.15) 

i=3:  ^03, +F(1)  +  F(-1)  =  2-(f2-F0)+^-«S3  (3.16) 

where  the  flux  is  represented  as 

F  =  F0  -v "  +  FI  -v]  +  F2  -v,^  +  F3  -vf  (3.17) 
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and  the  source  term  is  represented  as 


S  =  S0-v^^+S\-vl+  S2  -v^  +  S3  -vl 


Eqn.  (3.8)  can  be  rewritten  generally  as 


^  =  RHS 
dt 


(3.18) 


(3.19) 


An  explicit  time  update  of  eqn.  (3.19)  is  done  using  the  4*  order  Runge-Kutta  method.  The 
following  equations  show  how  the  time  update  is  done. 


e' ,  =a,+y'^«ia 


n+— 
2 


At 


Q  ,=Q„+--RHS 

«+2  2 


r  \ 

0*  1 

,  n+-  , 

V  2J 


QZ=Q„+At-RHS 


f  \ 

r, 

,  n+—  , 

V  IJ 


(3.20) 

(3.21) 


(3.22) 


0«+i  -  +  — 


r 

f  \ 

f  \ 

\ 

RHs{q)]^1RHS 

Q*  1 

+  2  •  RHS 

Q\ 

+rhs(q::) 

V 

n-\ — 

^  2  J 

n-\ — 

^  2j 

J 

(3.24) 


4.0  1-D  and  Quasi  1-D  test  cases 

Formulations  for  ID  Euler,  quasi-lD  Euler,  and  ID  magneto-gas  dynamics  (MGD)  equations  are 
presented  in  this  section. 

4.1  Perfect  Gas  Euler  Equations  in  1-D 

The  partial  differential  equations  (PDE)  that  result  from  the  conservation  principle  for  mass, 
momentum,  and  energy  are 
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Qt  Q- 


91 

P 

~fl 

pu 

2 

92 

= 

pu 

F  = 

f2 

p  +  pu 

.93  _ 

e 

L/3J 

(e  -h  p)u 

92 

2 

97 

P  +  — 

91 

(93  +  P)— 

91 


(4.1) 


where  /?  =  (y  -  l/e  -  —  pw 


2^ 
^  y 


=  (y  -1) 


93 


2  ^ 

l2_ 

291 


where  p  is  density,  u,  the  velocity,  p,  the  pressure,  e,  the  total  energy  per  unit  volume,  Q ,  the 
vector  of  dependent  variables  (conserved  quantities),  and  F  ,  the  flux  vector. 

4.2  Quasi  1-D  Perfect  Gas  Flow 

The  governing  equations  for  an  inviscid  quasi- ID  flow  are  given  by 

Qt+Fx=S 


Q  = 


■  puA  ■ 
/~  ~~2  1 . 

1 

_ 1 

to;'?! 

1 _ 

F  = 

\p  +  pu  lA 
(e  +  p^A 

S  = 

1 

1 _ 

Let  p  =  pA,  u=u,  e  =  eA,  p  =  pA,  then 


(4.2) 


"P" 

)  ‘’“2] 

0 

pu 

e 

F  = 

W  +  P"  1 

\e  +  pju 

S  = 

^2 

y3_ 

P-g 

0 

,  1  dA 

where  e  = - 

A  dx 
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where  A  is  the  area  of  the  nozzle.  That  is,  the  governing  equations  for  a  quasi- ID  flow  (eqn.  4.2) 
are  same  as  the  ID  Euler  equations  (eqn.  4.1)  except  for  the  fact  that  the  momentum  equation 
includes  a  source  term  and  the  fact  that  p,  e,  and  p  are  scaled  by  the  area  A. 

The  presence  of  the  source  term  S  necessitates  evaluation  of  the  following  integrals  in  the  weak 
form  of  the  governing  equations: 


h 

2  . 

\vg  ■  S2  •  dx ,  \  =  0,n-\  (4.3) 

h 

2 


This  integral  may  be  transformed  to 


|jv'  -s^-d^  ,  i  =  0,n-l 


(4.4) 


where  the  source  term  is  expanded  in  terms  of  the  basis  functions  as 

*2fe)=  S  . (4.5) 

2=0 

This  reduces  the  evaluation  of  the  source  term  integral  on  the  RHS  of  eqn.  (4.2)  to  the  evaluation 
of  the  integrals 


i=0.n.j=0.n 


The  final  form  of  the  source  term  in  eqn.  (3.18)  is 


SJ  = 


vOy 


(4.6) 


(4.7) 


4.3  Governing  Equations  for  MGD  Flows  in  1-D 

Simulation  of  inviscid  flow  of  a  perfect-gas  in  the  presence  of  a  magnetic  field  (inviscid  MGD 
flow)  involves  8  dependent  variables,  namely,  density  p,  velocity  components  (u,v,w),  energy 
per  unit  mass,  e,  and  three  components  of  the  magnetic  field,  (Bx,  By,  Bz).  In  the  case  of  one¬ 
dimensional  flows,  the  derivatives  in  two  coordinate  directions,  for  examples,  the  y  and  z 
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directions,  vanish  identically.  When  there  is  no  magnetic  field  present,  this  implies  that  the  v 
and  w  components  of  velocity  must  remain  constant  in  the  x-direction.  On  the  other  hand,  when 
a  magnetic  field  is  applied  to  the  flow  of  a  conducting  fluid,  Lorentz  forces  are  produced,  that 
can  cause  an  x-variation  of  v  and  w  components.  The  only  dependent  variable  that  vanishes 
identically  is  the  x  component  of  the  magnetic  field,  Bx.  (  The  equation  set  for  ideal  MHD  is 
over-determined,  and  requires  that  the  magnetic  field  is  divergence  free.  In  1-D  this  implies  that 

dB  dB^  ^  n  ^ 

=  constant ,  since  — ^  ^  =  0^  — ^  =  0 .  ) 

dy  dz  dx 

Further  details  on  the  conservation  laws  in  MOD  may  be  found  in  any  standard  text  on  the 
subject,  e.g.,  Sutton  and  Sherman  [5].  The  governing  equations  for  inviscid,  ID,  MGD  flows 
may  be  written  in  conservation  form  as 

Q,+Fx=0  (4.7) 

where  the  vector  of  dependent  variables,  Q,  and  the  flux  vector,  F,  are  given  by 


"  p  ’ 

pu 

pu 

pu^  -1-  p* 

pv 

puv-B^By 

pw 

F  = 

puw-B^B^ 

By 

ByU-B^v 

Bpi  -  B^w 

e 

[e  +  p^\i-B^  (b^u  +  ByV  +  B^w)_ 

Total  energy,  e,  and  modified  pressure,  p  are  defined  by 

c  =  ^  +  i  p(a’  +  +  w’ )+ 1(5/  +  S/  +  s;j 

P"  =P  +  \{Bx^+By^+By) 

Eigen  values  of  the  1-D  inviscid  flux  may  be  computed  from  the  above  relations,  and  are 
required  to  estimate  numerical  dissipation  in  the  Lax-Friedrichs  method.  They  are: 

“-‘■/I 

if  *2  .  2,  2 

—  a  +^Ja  -4a 
V  J 


X  = 


u  +  c^ 
u 


,  with 


(4.9) 

(4.10) 


(4.11) 
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b.  = 


*  = 


b.  = 


B. 


IrP 


b  =  Jb.^  +bj  +b.^ 


\v  »  r~2  Ti 

a  =  —  a  =  yj  a  +  b 

V  P 


4.4  Limiters 

In  the  numerical  solution  of  partial  differential  equations  that  govern  the  physics  of  fluid  motion 
limiters  are  a  necessary  evil.  While  they  do  reduce  the  accuracy  of  a  numerical  scheme,  they  are 
required  to  ensure  stability  when  discontinuities  and  high  curvature  regions  are  encountered. 
Even  though  Atkins  and  Shu  state  in  one  of  their  papers  [6]  that  they  have  been  able  to  obtain 
stable  solutions  for  a  non-linear  problem  in  ID  with  shocks,  their  solutions  for  DG-p2  exhibit 
typical  behavior  of  a  solution  obtained  without  limiters  (overshoots  and  undershoots  in  the  shock 
region)  clearly  confirming  the  need  for  limiters. 

Limiters  enter  into  a  computational  process  that  employs  a  piecewise  continuous  approximation 
for  the  dependent  variables  when  integrals  at  element  boundaries  are  evaluated.  Evaluation  of 
such  an  integral  involves  computation  of  the  integrant  at  element  boundaries  by  interpolating  the 
piecewise  continuous  approximation.  This  process  invariably  leads  to  two  different  values  for 
the  integrant  computed  from  the  interpolation  performed  for  the  two  neighboring  elements.  Such 
a  discontinuity  at  an  element  boundary  is  usually  resolved  using  a  Riemarm  solver.  Limiters  are 
required  to  ensure  that  the  interpolation  process  does  not  result  in  new  maxima  or  minima.  In 

our  formulation,  we  employ  a  minmod  limiter  [7]  for  each  of  the  coefficients  Qg{t)  in  eqn. 

(2.2).  More  research  is  needed  to  optimize  limiting  to  ensure  stability  of  the  scheme  without 
compromising  accuracy. 


5.0  1-D  Results 

Results  from  the  application  of  the  ID  DG  scheme  for  the  simulation  of 

1 .  a  simple  scalar  wave  propagation, 

2.  inviscid  propagating  density  wave, 

3.  inviscid  flow  in  a  shock  tube, 

4.  inviscid,  quasi- ID  flow  in  a  convergent-divergent  nozzle, 

5.  MGD  flow  with  an  intermediate  shock,  and 

6.  Brio-Wu  problem 

are  presented  in  the  following  sections.  Note  that  the  order  in  which  the  results  are  presented 
represent  the  increasing  order  of  difficulty  involved  in  simulating  a  simple  scalar  wave  to  an 
MGD  flow  that  involves  discontinuities. 

The  first  example  considered  here,  namely  propagation  of  a  simple  wave,  is  probably  the 
simplest  form  of  conservation  law  that  one  may  encounter.  In  this  case  the  vector  of  dependent 
variable  has  just  one  element.  In  other  words,  the  governing  equation  involves  a  scalar 
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conservation  law.  The  flux  vector  is  just  a  linear  funetion  of  the  dependent  variable  vector.  The 
solution  involves  a  simple  wave  propagating  at  a  eonstant  speed. 

The  second  example  involves  a  system  of  eonservation  laws.  The  equations  solved  are  the  one¬ 
dimensional  Euler  equations  (eqn.  4.1).  This  case  also  involves  a  simple  traveling  wave.  By 
setting  u=\  and  p  =  \,  the  solution  is  forced  to  eorrespond  to  a  density  wave  traveling  at  a 
constant  speed.  This  example  serves  the  purpose  of  validating  the  implementation  of  a  ID,  DG 
solver  for  Euler  equations. 

The  flow  in  a  shock  tube  is  a  elassieal  problem  employed  in  computational  fluid  dynamics 
(CFD)  to  evaluate  the  accuraey  of  a  numerical  scheme  in  simulating  a  system  of  equations  that 
involve  multiple  waves  traveling  at  non-eonstant  speeds. 

Inviscid,  quasi- ID  flow  adds  the  eomplexity  of  a  source  term  to  the  governing  equations  and 
serves  the  purpose  of  evaluating  the  aecuracy  of  a  numerical  scheme  in  dealing  with  souree 
terms. 

Intermediate  MHD  shoeks  are  diseontinuities  aeross  whieh  the  tangential  magnetic  field 
component  ehanges  sign  and  the  shock  frame  veloeity  of  the  flow  changes  from  super  to  sub- 
Alfvenic  across  the  shoek.  We  seleet  this  problem  for  study,  sinee  it  stands  to  benefit  from 
higher  order  simulations  and  possesses  smooth  solutions  (when  used  with  periodie  BCs  in  1-D) 
until  the  formation  of  the  shoek. 

The  Brio  and  Wu  [2]  1-D  MHD  shoek  tube  validation  ease  represents  an  important  porting  of 
shoek  tube  validation  studies  (sueh  as  that  of  Sod,  Van  Leer,  etc.  in  perfect  gases,)  to  MHD. 


5.1  Scalar  Wave  Propagation  in  ID 


This  simple  time-dependent  problem  demonstrates  the  hierarchieal  improvement  in  simulation 
accuraey  attainable  using  the  DG-pO,  DG-pl,  and  DG-p2  sehemes.  In  this  ease  the  governing 
equation  is  given  by 

-0.1<x<2.0  (5.1) 


with  a  cubie  polynomial  as  the  initial  condition.  That  is, 


0, 


u{x,0)  = 


- 16000 +1200-x^  +1 
- 16000 -x^  +1200-x^  +1 


-0.1  <x< -0.05 


-0.05  <x<0 


0  <  X  <  0.05 


0 


0.05<x<2.0 


(5.2) 
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The  boundary  conditions  employed  are 


w(-0.1,t)  =M;(;(-0.1,t)  =u{2,t)  =  =  0 


(5.2) 


Cubic  Pulse  (10  points) 


Figure  1.  Solution  for  t=  1.6 


The  grid  employed  has  a  resolution  of  dx  =  0.01.  Solutions  obtained  for  t  =  1.6  using  DG-pO, 
DG-pl,  and  DG-p2  are  compared  with  the  exact  solution  in  fig.l.  The  exact  solution  is  nothing 
but  a  simple  wave  traveling  to  the  right  at  a  constant  speed.  As  such  the  solution  at  t  =  1.6  is  the 
same  as  the  solution  at  t  =  0  shifted  to  the  right  by  1.6 . 

Hierarchical  nature  of  the  simulation  accuracy  may  be  seen  from  the  fact  that  while  the  DG-p2 
solution  agrees  with  the  exact  solution  almost  exactly,  DG-pl  solution  is  visibly  different  from 
the  exact  solution  and  DG-pO  solution  has  an  unacceptably  large  error.  The  dissipative  and 
disspersive  nature  of  the  errors  may  be  observed  from  the  decrease  in  amplitude  and  increase  in 
wave  length  that  accompany  a  decrease  in  the  order  of  the  scheme. 

5.2  Propagating  Density  Wave 

In  the  last  example,  we  demonstrated  the  accuracy  of  a  DG  scheme  in  simulating  a  scalar 
traveling  wave.  The  governing  equation  involved  just  one  PDE  (partial  differential  equation). 
Here  we  consider  a  system  of  PDFs  representing  three  conservation  laws,  namely,  the  ID,  Euler 
equations  (eqn.  4.1).  In  order  to  simplify  the  problem  and  employ  this  case  as  a  test  case  for  the 
verification  of  our  implementation  of  a  DG  scheme  for  ID,  Euler  equations,  we  force  u  and  p  to 
be  identically  equal  to  1  so  that  the  solution  involves  a  simple  density  wave  propagating  at  a 
constant  speed.  In  reality  this  case  is  no  different  from  the  first  example  as  far  as  the  solution  is 
concerned.  The  only  difference  is  that  we  use  a  code  developed  for  solving  the  ID,  Euler 
equations. 

The  ID,  Euler  equations  are  solved  in  the  domain  -l<x<l,  with  the  initial  condition  for 
density  given  by 
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p  (x,0)  =  2  +  sin(7u:) 

With  periodic  boundary  conditions,  u  =1,  and  p  =  l  this  problem  has  the  following  exact 
solution: 


p  =2  +  sin(7i  (x  -  Mt)) 

The  solution  for  this  case  at  t=100  is  shown  in  fig.  2.  It  may  be  seen  from  this  figure  that  the  p2 
scheme  preserves  the  wave  better  than  the  pi  scheme  indicating  the  improved  accuracy  that  is 
possible  with  the  p2  scheme. 

Error  analysis  performed  for  three  different  grid  spacings  is  shown  in  Table  1.  It  may  be  seen 
from  the  table  that  the  error  for  the  p2  scheme  is  2  to  3  orders  of  magnitude  smaller.  The  order 
of  convergence  is  also  higher  for  the  p2  scheme. 


No.  of  grid  points 

PI:  error 

order 

P2:  error 

order 

32 

2.5e-l 

2.0e-3 

64 

3.5e-2 

2.8 

l.le-4 

4.1 

128 

4.5e-3 

2.9 

9.6e-6 

3.6 

Table  1:  Error  for  the  pi  and  p2  schemes  for  different  grid  resolutions 


^ine  wavi0 


Figure  2.  Solution  for  a  simple  density  wave  at  t=100. 

One  of  the  major  benefits  that  one  may  expect  to  derive  from  application  of  a  higher  order 
scheme  to  an  unsteady  problem  involves  predicting  solutions  for  large  t.  In  order  to  assess  the 
improvement  in  accuracy  attainable  with  a  fourth  order  scheme  (DG-p3),  the  density  wave  was 
propagated  to  t= 5 0,000.  The  results  obtained  using  the  p2  and  p3  schemes  are  compared  in 
fig.3. 
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The  grid  resolution  employed  corresponds  to  20  points  per  wave.  It  may  be  seen  from  this  figure 
that  the  p3  scheme  performs  much  better  than  p2  scheme  justifying  the  use  of  higher-order 
accuracy  in  similar  problems. 


Density  Sine  Wove 


)  ^  ■ _ I _ I _ i _ I _ I _ _ I _ i 

0.0  0.2  0.4  0.6  0.8  1.0 


X 

Figure  3.  DG-p2  and  DG-p3  schemes  compared. 


5.3  Shock  Tube 

The  flow  in  this  example  starts  with  a  stationary  gas  in  a  tube  partitioned  by  a  diaphragm.  The 
gas  on  the  right-hand  side  of  the  diaphragm  is  at  a  higher  pressure  {pA)  and  density  (p4)  than 
the  gas  on  the  left-hand  side  {p\  and  pi) .  Simulations  were  carried  out  with  pressure  and  density 
ratios  given  by 


p\  pi 

When  the  diaphragm  is  punctured  at  t  =  0 ,  a  shock  wave  propagates  into  the  high  pressure 
region  while  an  expansion  wave  and  a  contact  discontinuity  propagate  into  the  low  pressure 
region.  The  simulation  discussed  here  employed  a  computational  domain  given  by: 
-  0.5  <  X  <  0.5 .  The  diaphragm  was  located  at  x  =  0 .  The  grid  employed  had  300  elements. 

Solutions  obtained  using  DG-pO,  DG-pl,  and  DG-p2  are  compared  with  the  exact  solution  in 
fig.4.  While  the  figure  shows  a  substantial  improvement  in  accuracy  as  the  order  of  the  scheme 
is  increased  from  1  (DG-pO)  to  2  (DG-pl),  the  solution  for  3^“^  order  (DG-p2)  scheme  appears  to 
be  no  better  than  the  solution  for  the  2"‘*  order  scheme.  This  may  be  due  to  the  fact  that  unlike 
the  last  two  examples  the  2“*^  derivative  in  this  case  is  identically  zero  for  a  substantial  portion 
of  the  computational  domain  decreasing  the  role  of  the  parabolic  term  in  the  basis  function. 
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Fig.  5  shows  the  effect  of  the  CFL  number  on  the  solution  for  the  DG-pl  scheme.  While 
CFL=0.1  and  CFL=0.3  produce  almost  identical  solutions,  CFL=0.4  results  in  an  oscillatory 
behavior  indicating  that  the  simulation  is  unstable.  This  behavior  is  consistent  with  the  fact  that 
the  CFL  limit  for  a  DG-pn  scheme  is  given  by  [8] 


CFL< - 

2-n  +  \ 


(37) 


and  as  such  the  DG-pl  scheme  is  unstable  for  CFL  >  — . 

3 


Shock  Tube  P4/P1=2,  RH04/RH01=2  Shock  Tube  P4/P1=2.  RH04/RH01=2 

Oisconlinuous  Goterkin  1 D  Euler  Code  Discontinuous  Golerkin  1 D  Euler  Code 


Figure  4.  DG-pO,  DG-pl,  and  Dg-p2  compared  Figure  5.  DG-pl;  effect  of  CFL 


5.4  Quasi-ID  Flow 

So  far  we  have  considered  a  simple  wave  resulting  from  a  scalar  conservation  law,  a  density 
wave  that  represents  the  solution  to  a  system  of  conservation  laws,  and  a  flow  that  involves 
multiple  waves  traveling  at  non-constant  speeds.  In  this  section  we  present  solutions  for  inviscid 
flow  in  a  quasi- ID  convergent-divergent  nozzle  which  involves  solving  a  system  of  conservation 
laws  with  an  added  complexity  of  a  source  term. 

The  convergent-divergent  nozzle  selected  for  this  example  has  cross-sectional  area  represented 
by  a  cubic  polynomial.  The  computational  domain  extends  from  x  =  -l  to  x  =  l.  The  area 
distribution  is  symmetric  about  x  =  0 .  Area  is  minimum  at  x  =  0 .  The  area  ratio,  that  is  the  ratio 
between  the  area  at  x  =  0  and  area  at  x  =  -1  (or  x  =  1)  is  equal  to  0.9.  The  nozzle  contour  was 
defined  in  such  a  way  that  the  slope  was  zero  at  the  entrance  (x  =  -1),  exit  (x  =  1),  and  throat 
(x  =  0)  of  the  nozzle.  Flow  conditions  were  chosen  to  ensure  that  the  flow  did  not  choke,  that  is 
the  flow  remained  either  supersonic  or  subsonic  everywhere  depending  on  the  flow  at  nozzle 
entrance. 
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The  results  obtained  from  our  study  are  shown  in  fig.6.  The  DG-pO  scheme  has  a  slope  of  1 
while  DG-pl  has  a  slope  of  1.92  which  corresponds  to  orders  of  accuracy  of  1  and  1.92 
respectively.  We  have  not  yet  been  able  to  obtain  satisfactory  performance  with  the  DG-p2 
scheme.  The  error  for  the  DG-p2  scheme  (not  shown  in  the  figure)  is  slightly  less  that  for  DG-pl 
but  only  slightly.  We  believe  that  this  is  due  to  the  fact  that  the  error  for  the  DG-pl  solution  is 
already  very  near  machine  zero  for. 


Supersonic  CD  Nozzle  AMIN/AMAX=0.g  M=2.07 


Figure  6.  Error  Analysis 

5.5  Intermediate  Shock  Simulation 

Intermediate  MHD  shocks  are  discontinuities  across  which  the  tangential  magnetic  field 
component  changes  sign  and  the  shock  frame  velocity  of  the  flow  changes  from  super  to  sub- 
Alfvenic  across  the  shock.  Wu  [3]  first  demonstrated  the  feasibility  of  formation  of  intermediate 
shocks  in  a  MHD  flow  through  nonlinear  wave  steepening  from  continuous  waves.  We  select 
this  problem  for  study,  since  it  stands  to  benefit  from  higher  order  simulations  and  possesses 
smooth  solutions  (when  used  with  periodic  BCs  in  1-D)  until  the  formation  of  the  shock.  Wu  [3] 
employed  the  following  simple  wave  solution  of  the  MHD  equations  as  the  initial  conditions: 


=  0.5  •  sin(27u:) 
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dp  ^ 

Cj-a^ 

w  =  0  =  I  =0 

J 

Jiang  and  Wu  [9]  obtained  numerical  solutions  for  this  intermediate  shock  problem  using  a  3 
order  ENO  and  a  5*  order  WENO  schemes.  Their  results  showed  that  the  3’^'*  order  ENO  and  the 
5*  order  WENO  schemes  performed  quite  well  when  marched  in  time  to  obtain  solution  for 
t=0.25  (fig.  7).  At  this  particular  instance  in  time  the  y-component  of  the  magnetic  field,  By,  is 
quite  smooth  and  still  far  from  steepening.  On  the  other  hand,  at  t=0.5  (fig.  8),  By  profile 
exhibits  large  gradients  showing  a  tendency  to  form  a  shock.  Their  error  analysis,  when  applied 
to  the  whole  computational  domain  showed  poor  results.  When  they  restricted  their  error 
analysis  to  smooth  regions  of  the  flow,  they  were  able  to  show  that  their  algorithm  had  the 
desired  order  of  accuracy. 

Solutions  obtained  from  our  numerical  simulation  for  this  problem  are  shown  in  figs.  7-8. 
Solutions  have  been  obtained  using  various  grids  with  spatial  resolution  as  coarse  as  0.1  (10 
elements)  to  as  fine  as  0.00039  (2560  elements).  The  y-component  of  the  magnetic  field  at 
t=0.25,  shown  in  fig.  7,  appears  to  have  been  quite  well  predicted  by  all  grids  except  the 
coarsest.  The  plot  of  the  LI -norm  of  the  error,  computed  using  the  finest  grid  solution  as  the 
exact  solution,  indicates  that  the  2“^*  (pi)  and  3'^'^  order  (p2)  schemes  have  the  same  level  of 
accuracy  and  match  the  accuracy  of  Jiang  and  Wu’s  order  ENO  scheme.  The  error  analysis 
for  t=0.5  indicates  that  the  3’^'*  order  (p2)  DG  algorithm  is  more  accurate  than  the  2"*^  order 
scheme.  More  interestingly,  both  the  2“^*  and  3"^^*  order  DG  algorithms  perform  better  than  even 
the  5*  order  WENO  scheme  [9]. 

The  fact  that  the  p2  solution  is  more  accurate  than  pi  solution  for  t=0.5  and  the  fact  that  for 
t=0.25  the  two  solutions  have  almost  the  same  order  of  accuracy  is  probably  due  to  the  fact  that 

the  contribution  from  Vg(^)  is  not  significant  when  the  solution  is  smooth  and  becomes 
significant  as  the  By  profile  starts  steepening. 


Figure  7.  Solution  for  t=0.25 
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Figure  8.  Solution  for  t=0.5 

5.6  The  Brio-Wu  Shock-Tube  Problem 

Brio  and  Wu  [2]  presented  the  first  implementation  of  Roe’s  Riemann  solver  in  1-D  MHD 
problems  and  obtained  a  Roe  linearization  for  the  case  of  y  =  2.  The  shock  tube  problem  used  by 
them  to  demonstrate  their  solver  represents  an  important  porting  of  shock  tube  validation  studies 
(such  as  that  of  Sod,  Van  Leer,  etc.  in  perfect  gases,)  to  MHD. 


p=  1  1 

p=-i 

P=1  1 

p  =  0.125 

By=-1 

o 

II 

N 

PQ 

\  Bz=  0 

u,v,w  =  0  1 

\  u,v,w  =  0 

Diaphragm 


Figure  9.  Initial  Conditions  for  the  shock- tube  problem 


The  initial  conditions  employed  in  this  problem  are  shown  in  fig.  9.  Distributions  of  pressure, 
density,  and  y-component  of  the  magnetic  field  obtained  from  our  numerical  simulation  using  a 
fine  grid  (800  elements)  and  a  coarse  grid  (100  elements)  are  shown  in  figs.  10-11.  Our  fine  grid 
solutions  using  2“*^  order  (pi)  and  order  (p2)  schemes  were  compared  with  those  from  a 
Kinetic  Flux  Vector  Splitting  scheme  developed  by  Xu  [10]  and  implemented  earlier  in  MHD 
related  work  in  HyPerComp  (Munipalli  and  Shankar  [11]).  These  two  results  match  closely. 
Coarse-grid  solutions  are  compared  with  the  3"^^*  order  fine-grid  solution  in  Fig.  11.  It  may  be 
seen  from  this  figure  that  the  2“^*  and  order  solutions  are  almost  identical  and  show  a  marked 
improvement  over  the  L*  order  solution.  As  in  the  case  of  intermediate  shock  solution  for 
t=0.25,  the  2“'^  derivatives  are  almost  non  existent  in  many  regions  of  the  flow  contributing  to  the 
fact  that  the  2“*^  and  order  solutions  are  almost  identical. 
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Figure  10.  Fine-grid  (800  grid  elements)  solution  for  t=0.4. 


Figure  11.  Coarse-grid  (100  grid  elements)  solution  for  t=0.4. 


6.0  Point  Implicit  Scheme  with  Orthogonal  Polynomials 

All  the  computations  discussed  so  far  were  carried  out  using  a  4*  order  Runge-Kutta  explicit 
time  marching  scheme.  The  need  for  an  implicit  scheme  arises  from  the  fact  that  the  CFL  limit 
for  an  explicit  DG-pn  scheme  is  given  by  [8] 


CFL<—^  (6.1) 

ln  +  \ 

This  imposes  a  very  stringent  limitation  on  the  allowable  time-step.  Implicit  schemes,  on  the 
other  hand,  do  not  have  any  CFL  limit  and  the  allowable  time  step  is  determined  only  by  the  time 
accuracy  required  by  the  problem  under  consideration.  For  steady  state  simulations,  ability  to 
employ  large  values  for  the  time-step  will  result  in  fast  tum-around  time.  Also,  when  a  non- 
uniform  grid  is  employed,  allowable  time-step  is  determined  by  the  finest  grid  element  which 
could  be  orders  of  magnitude  smaller  than  the  coarsest,  especially  for  viscous  flows.  Hence,  it  is 
highly  desirable  to  develop  an  implicit  scheme. 

In  order  to  advance  the  solution  from  time  level  n  to  level  {n+1)  using  an  implicit  scheme,  all  the 
terms  in  the  discretized  governing  equations  (eqn.  3.8)  must  be  evaluated  at  time  level  {n+1). 
Specifically,  all  the  flux  terms  and  the  source  terms  need  to  be  computed  at  time  level  {n+1). 
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Expressions  for  F”’'’^(i)=  F 


andF"+^(-l)=F 


^-h  ^ 


V  ^  y 


in  eqn.  3.8  (flux  at  the  right 


V  ^  y 


and  left  boundary  respectively)  using  the  scalar  dissipation  scheme  (eqn.  2.6)  and  orthogonal 
polynomials  as  basis  functions  is  presented  in  this  section. 

6.1  Derivation  of  the  scheme 
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PDE  for  QO : 
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dt  I  2  2  j 


From  equations  (6.2),  (6.3),  and  (6.4),  we  get 
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From  equations  (6.2),  (6.3),  and  (6.4),  we  get 


he  dQlc  ^JheY  h^r 

2  dt  2  2 


+  Ato  - 


+  At|  (^1  ^  |  +  a|  ^  + At|  ,1,1  ^  |-a| 

^  2  dt 


hr  Y  dQOfir  dQlrir 


hc\\5Q^c 

2]]  dt 


+  +  ^  U2.FO 

I  2  j  dt  dt 


h.c  =  —\F\^-F\-^-2-FQ 

hA  [2]  [2] 


24 


(6.7) 
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We  employ  a  point  implicit  scheme  to  solve  for 


dt 


dQ\ 

dt 


from  equations  6.6  and  6.7. 


Point-implicit  scheme  is  an  iterative  scheme  in  which  the  value  from  the  previous  iteration  is 
used  for  all  the  terms  in  an  equation  except  the  term  that  corresponds  to  the  dependent  variable 
that  is  solved  for.  The  convergence  characteristics  of  such  a  scheme  for  numerical  simulation  of 
a  quasi- ID  flow  in  a  convergent  divergent  nozzle  is  compared  with  the  convergence 
characteristics  of  a  four  stage  explicit  Runge-Kutta  time-marching  scheme  in  fig.  12.  The  figure 
shows  that  an  order  of  magnitude  acceleration  in  convergence  may  be  achieved  using  a  point 
implicit  scheme. 


Figure  12.  Convergence  Characteristics  of  a  pi  Point  Implicit  Scheme  for  Quasi- ID  Flow  in  a 

Convergent-Divergent  Nozzle 

7.0  DG  Formulation  for  the  3D  Euler  Equations 

In  this  section,  the  formulation  of  the  3D  DG  Euler  equations  is  presented.  The  structure  of  the 
equations  follows  the  notation  and  conventions  used  in  the  3D  CEM  DG  code  TEMPUS. 

The  3D  Euler  equations  can  be  written  as 
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(7.1) 


(7.2) 


where  p  is  density,  u  is  the  velocity  in  the  x-direction,  v  is  the  velocity  in  the  y-direction,  w  is 
the  velocity  in  the  z-direction,  e  is  the  total  energy  per  unit  volume,  and  p  is  the  pressure  defined 
as 


/?  =  (y  -l{e-^p{u^  +  where  y  is  the  ratio  of  specific  heats. 


The  DG  form  of  the  Euler  equations  can  be  written  as 


|v '  =  -  |v  ‘f(q)- ms  +  J Vv  '  •  f[^V 


(7.3) 


where  ic  is  the  current  cell,  v'  the  i*  basis  function,  V  the  cell  volume,  and  S  the  cell  face 
surface. 


The  basis  functions  in  this  version  of  the  code  are  monomials.  The  monomials  for  up  to  P=2 
(third  order)  are  given  in  eqn.  (7.4).  For  first  order,  P=0,  i=0,  and  for  second  order,  P=l, 
0<i<3. 

^0=1 

Vi=x-x^,V2=y-y„V3=z-z^ 

V4  =  (x - xj- (y - yj,  V5  =  (x - xj- (z - zj,  v,  =  (y - yj- (z - zj 
V7  =  (^ - Vs  =  (t  - Tc)".  V9  =  (z  - zJ 

where  c  is  for  the  cell  centroid. 


The  solution  vector  Q  is  expanded  in  the  basis  functions  and  the  component)  of  the  vector  Q  is 
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(7.5) 


JP 


1=0 


with  JP=0  for  P=0,  JP=3  for  P=l,  and  JP=9  for  P=2. 


Similarly,  the  flux  vector  F  is  expanded  in  a  Taylor  series  about  the  cell  centroid.  The  expansion 
for  the  component  of  the  flux  vector  may  be  written  as 

F,=fv'f‘  (7.6) 

/=0 


with 


f;=Fj(q^)  for  i=0  (P=0) 


(7.7) 


5  QP 

fori=l,3  (P=l) 

k=l 

f] = 

for  i=4,  il=l,  i2=2,  for  i=5,  il=l,  i2=3,  for  i=6,  il=2,  i2=3,  for  i=7,  il=l,  i2=l,  for  i=8,  il=2, 
i2=2,  and  for  i=9,  il=3,  i2=3. 

8.0  Evaluation  of  Integrals  in  eqn.  7.3 

This  section  describes  how  each  of  the  integral  terms  in  eqn.  (7.3)  is  evaluated.  The  first  integral 

to  be  evaluated  is  [v '  -^-dV  which  will  be  rewritten  using  the  basis  function  expansion  for  Q. 

'  dt 
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d^dr\di^  ,  [v '  -^-dV  can  be  rewritten  in  transformed  coordinates  as 
•’  dt 
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M  =  \my  \  is  the  mass  matrix. 
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The  next  integral  to  be  evaluated  is  j  Vv '  •  . 


W 


Following  the  evaluation  process  of  |v '  ^^-dV ,  J  Vv '  •  fi^^dV  is  rewritten  using  the  basis 
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Finally  using  the  expression  for  dV  given  earlier,  we  get 
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j  Vv '  •  F[Q)dV  =  ^(mxyfj  +  my^jfj  +  mz-jf^  (8.2) 

V-  7=0 

where  Mx,  My,Mz  =  |  |  J  ^^e  the  gradient  mass  matrices. 

The  final  integral  is  Jv  'f{^-  hdS  .  This  is  a  surface  integral  at  the  interior  cell  faces  which  will 

Sic 

use  information  from  the  current  cell  (ic)  and  the  neighbor  cell  (inb).  The  face  flux  is  formed 
using  the  local  Lax-Friedrichs  scheme. 

^(e)- " = I (4 + )+ ^  V..  (a  -  a* ) 

where  =  Iw  +  cl ,  u  =  u  ■  n+v  ■  n^,  +  w  -  n  and  c  = 
n  =  (n^,ny,n^'^  is  the  unit  normal  at  the  cell  face  and 
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F.^  =  fj,-n^  +  fy-hy+f^-h^  (all  expansions  from  cell  centroid  ic) 

^inb  -  fx'^x'^  fy'^y'^  fz'^z  cxpaosions  from  cell  centroid  inb) 

Expand  F  and  Q  in  the  basis  functions  and  jv  'f{^- MS  can  be  rewritten  using  the  following 
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expression  for  dS  = 
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Finally,  ^v'f{^-MS  can  be  rewritten  as 


where  S  =  [5,^.  J  is  the  surface  matrix  and  F  =  \fy  \  is  the  face  matrix. 

At  boundaries,  the  flux  computation  follows  the  same  structure  except  that  now 

Qinb  -  fiinction  )  with  the  functional  dependency  determined  by  the  boundary  condition 

employed.  Currently  inviscid  surface,  supersonic  inflow,  and  a  combination  of  subsonic 
inflow/outflow  boundary  conditions  have  been  implemented. 

Now  that  all  the  integrals  have  been  evaluated,  the  3D  DG  Euler  equations,  eqn.  (7.3),  can  be 
written  as 


=  +m«)+(*fc'Z  +Myl  +Mz-l\  (8.4) 
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All  the  integrals  in  M,  S,  F,  Mx,  My,  and  Mz  are  evaluated  using  Gaussian  quadrature.  The 
solution  is  advanced  in  time  using  a  4  stage  Runge-Kutta  time  algorithm. 

9.0  Conversion  of  CEM  Solver  to  Euler  Flow  Solver 

The  structure  of  the  CEM  solver  is  presented  in  the  form  of  a  flowchart  in  Appendix  1 .  As  was 
stated  in  our  proposal,  it  has  been  possible  to  convert  the  CEM  solver  to  an  Euler  solver  without 
making  any  appreciable  change  to  the  structure  of  the  solver.  The  details  of  the  Euler  solver  are 
presented  in  the  following  paragraphs. 

Basis  function/Geometry  Data  Base 

The  subroutines  hex_yol_moms.c,  hex_sur_moms.c,  and  hex  J'acjnoms.c  are  used  in  original 
form  to  compute  the  M,  Mx,  My,  Mz,  S,  and  F  matrices  which  will  be  used  in  solving  eqn  (8.4). 

Flow  initialization 

The  initialization  subroutine  qjnit.f  has  been  modified  to  initialize  the  flow  field  using  a 
constant  freestream  flow  condition  defined  by  the  user  in  the  form  of  the  variables  uincqn,  n 
=1,6.  This  initialization  needs  to  be  recoded  later  for  more  general  flow  initialization. 

Flux  computation 

The  flux  subroutine  iflxn.f  was  modified  to  compute  the  Euler  fluxes  using  eqn  (8.3).  This 
subroutine  is  modified  only  slightly  from  the  CEM  version.  In  CEM,  the  fluxes  fc  and  fnb  are 
linear  combinations  of  qic  and  qinb-  In  the  case  of  inviscid  flow,  eqns.  (7.6)  and  (7.7)  need  to  be 
used  since  the  fluxes  are  nonlinear.  Once  the  flux  coefficients  f-'  are  computed  using  eqns.  (7.6) 
and  (7.7),  the  original  structure  of  the  subroutine  can  be  retained  since  the  flux  subroutine  was 
developed  using  the  linear  combination  of  the  solution  coefficients  q^  as  the  flux.  The  S  and  F 
matrices  in  this  subroutine  are  used  without  change  from  the  original  subroutine.  Subroutines 
jacobian.f  and  fluxcoefficient.f  'wqvq  added  to  compute  the  flux  coefficients. 

Boundary  conditions 

In  this  version,  boundary  condition  subroutine  ubc _pc.f  has  been  modified  to  be  a  inviscid 
surface  tangency  boundary  condition,  ubc_rc.f  has  been  modified  to  be  a  plane  of  symmetry 
boundary  condition,  and  ubcjyb.fhas  been  modified  to  be  a  subsonic/supersonic  inflow/outflow 
freestream  boundary  condition.  New  subroutines  hex_ob_quads.c  and  hex_rc_quads.c  were 
added  (essentially  just  a  copy  of  hex _pc_quads.c)  with  the  appropriate  “oh”  or  “rc”  variables 
dynamically  allocated  and  used.  These  new  subroutines  allow  these  boundary  condition 
subroutines  to  use  the  same  data  base  that  the  “pc”  subroutine  is  using.  A  more  complete 
boundary  condition  library  needs  to  be  developed  since  a  general  flow  code  needs  many  more 
types  of  boundary  conditions. 

Gradient  Mass  Term 

The  subroutine  zv/wx./ needs  to  be  modifed  for  the  Euler  fluxes  (again  since  in  CEM  the  fluxes 
are  a  linear  combination  of  Q).  This  subroutine  computes  the  gradient  mass  integral  term  of  eqn. 
(8.2).  Again  where  solution  coefficients  were  being  used  for  the  flux,  flux  coefficients  are 
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computed  and  used.  This  is  just  the  same  procedure  as  was  done  in  modifying  the  flux 
computations. 


Post-processing 

Finally  the  output  subroutine  write _soln.f  to  be  modifed  to  output  the  appropriate  flow 

variables  for  the  Euler  flow  solver  and  remove  all  CEM  output. 

In  the  eonversion,  one  can  see  that  the  modularity  of  the  CEM  TEMPUS  code  allows  most  of  the 
code  to  be  used  in  the  original  state.  Only  four  segments  of  the  eode  needed  to  be  modified, 
namely,  (1)  flow  initialization;  (2)  flux  computation;  (3)  boundary  conditions;  and  (4)  post¬ 
processing. 

10.0  Verification  of  the  Implementation 

In  this  section,  we  present  results  obtained  using  the  3D  DG  Euler  solver  for  inviscid  flow  over  a 
2D  cylinder.  The  grid  employed  in  the  simulation  is  shown  in  fig.  13.  The  geometry  considered 
is  a  half  eylinder  with  radius  =1.0,  and  span  width  =  10.0.  The  outer  boundary  is  a  half  eirele  of 
radius  =  30.0.  The  grid  has  30  points  equally  spaced  on  the  circumferenee  of  the  cylinder  and  30 
points  normal  to  the  cylinder  surface  with  the  first  point  off  the  cylinder  surface  at  a  distance  = 
0.05.  Four  equally  spaced  span  points  are  used.  The  cell  topology  is  hexagonal.  This  grid  was 
run  using  the  3D  DG  Euler  eode  and  the  fully  validated  structured  grid  solver,  USA  [12].  The 
flow  conditions  employed  are 

p=1.0,  /7  =  1.0,  M  =  0.3,  V  =  0.0,  w  =  0.0,  y=1.4 


Figure  13.  2D  cylinder  grid 


The  DG  code  was  run  in  pO  (first  order)  mode  with  a  CFL=1.0  and  in  pi  (seeond  order)  mode 
with  CFL=0.3.  Fig.  14  shows  that  the  DG  code  is  slightly  less  dissipative  than  the  USA  code  for 
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the  same  order  of  aceuraey.  The  potential  flow  solution  is  shown  as  a  referenee.  One  must 
remember  that  in  solving  the  Euler  equations,  there  will  be  numerieal  dissipation  whieh  will 
cause  differences  from  the  potential  flow  solution.  Further  testing  of  the  Euler  DG  code  needs  to 
be  done.  Validation  of  the  p2  (third  order)  mode  will  be  carried  out  in  Phase  II. 


2D  Euler  Cylinder 


0.0  60.0  120.0  180.0 
Thelo  (Deg) 

Figure  14.  u  velocity  on  cylinder  surface 

11.0  Design  of  a  Common  Platform  for  Multi-physics  Computation 

In  this  section,  we  discuss  the  design  of  a  common  platform  for  multi-physics  computation. 

11.1  Multidisciplinary  analysis 

In  a  multidisciplinary  analysis,  from  the  point  of  view  of  numerical  simulation,  following 

situations  may  arise: 

1.  Different  physical  phenomena  involved  in  the  analysis  may  require  different 
computational  domains  (Fig.  15).  For  example,  in  the  case  of  aeroelasticity  the 
computational  domains  for  the  equations  of  fluid  flow  and  elasticity  are  different.  The 
two  phenomena  interact  at  the  boundaries  of  the  domain. 

2.  The  physical  phenomena  involved  may  require  the  same  computational  domain  (Fig.  16). 
In  this  case  there  are  two  different  possibilities: 

a.  The  same  grid  may  be  employed  in  solving  all  the  governing  equations. 

b.  All  physical  phenomena  involved  may  not  require  the  same  grid  resolution  and  as 
such  different  grids  may  be  used  for  different  physical  phenomena. 

3.  A  combination  of  situations  1  and  2. 
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Grid  for  Physics  1 
Grid  for  Physics  2 


Grid  for  Physics  3 
Grid  for  Physics  4 


Figure  15.  Different  computational  Figure  16.  Same  computational  domains 
domains  for  different  physics  for  different  physics 

In  all  these  situations  the  governing  equations  may  either  be  solved  in  a  coupled  manner  or  in  an 
uncoupled  manner.  A  coupled  solver  would  involve  inversion  of  a  single  matrix  while  an 
uncoupled  solver  would  involve  inverting  several  matrices.  In  either  case  a  highly  accurate 
interpolation  module  is  required  when  the  volume  grids  in  the  same  computational  domain  or 
surface  grids  in  the  case  of  different  computational  domains  are  different.  The  DG  scheme 
makes  such  interpolations  quite  simple.  The  expression  for  the  expansion  of  the  dependent 
variables  in  terms  of  the  basis  functions  also  serves  the  role  of  interpolation  function.  The 
TEMPUS  code  has  a  modular  structure  to  be  able  to  be  modified  to  handle  the  three  situations 
listed  above. 

The  partial  differential  equations  that  govern  different  physical  phenomena  may  either  be 
hyperbolic,  parabolic  or  elliptic.  From  a  numerical  simulation  point  of  view  they  may  be 
categorized  as  either  in  a  conservation  form  or  in  a  non-conservation  form.  In  the  case  of 
equations  in  conservation  form  with  all  positive  eigen-values,  the  solution  process  mainly 
involves  computation  of  flux  vectors  at  element  boundaries,  volume  source  terms,  and  in  the 
case  of  implicit  schemes,  the  flux  Jacobian.  Therefore,  it  is  indeed  possible,  as  has  been 
demonstrated  by  us  in  the  case  of  computational  electromagnetics  (CEM)  and  computational 
fluid  dynamics  (CFD),  that  a  common  tool  may  be  developed  with  individual  modules  for 
computation  of  flux  terms,  source  terms  and  flux  Jacobian  to  account  for  different  physical 
phenomena  involved. 

The  discussion  above  leads  to  the  following  ingredients  for  a  multidisciplinary  tool: 
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1.  A  module  for  establishing  eonneetivity  at  the  boundaries  of  different  computational 
domains  (such  as  solid/liquid  interface). 

2.  A  routine  to  compute  connectivity  between  different  grid  cells  in  the  same  computational 
domain  but  belonging  to  different  grids. 

3.  A  solver  for  equations  in  conservation  form  with  only  positive  eigen- values  that  includes 
individual  flux,  source,  and  Jacobian  modules  for  various  physical  phenomena  involved 
in  the  simulation. 

4.  Individual  solvers  for  all  other  governing  equations. 

5.  Interface  for  exchanging  information  at  computational  boundaries  for  different  physics. 

In  addition  to  having  the  ingredients  listed  above,  a  successful  commercial  multidisciplinary 
simulation  tool  should  be  complete,  user  friendly,  fast,  and  should  require  reasonable  computing 
resources.  Completeness  requires  the  tool  to  perform  all  the  tasks  related  to  the  simulation.  User 
friendliness  is  obtainable  through  the  development  of  a  graphical  user  interface  (GUI)  that  is 
intuitive  and  that  guides  the  user  through  the  process  quite  smoothly.  Speed  and  resource 
requirements  are  closely  associated  with  the  numerical  scheme  employed  in  the  simulation  and 
the  extent  of  intelligence  built  into  the  tool. 

We  will  use  HyPerComp’s  CEM  simulation  tool  TEMPUS  as  the  initial  framework  to  develop 
this  multi-disciplinary/multi-physics  tool.  The  3D  DG  CEM  solver  part  of  TEMPUS  has  been 
mentioned  previously  in  the  report  as  it  is  the  solver  which  was  used  in  the  conversion  to  a  3D 
DG  Euler  solver.  Now  the  whole  TEMPUS  environment  will  be  described. 


11.2  TEMPUS 

In  the  last  three  years,  HyPerComp  has  been  involved  in  the  development  of  TEMPUS,  a 
complete  tool  for  the  simulation  of  Maxwell’s  equations.  Starting  with  the  creation  of  the 
geometry  of  the  system  to  be  simulated,  TEMPUS  has  the  ability  to  repair  geometric 
imperfections,  define  a  computational  domain  and  discretize  it,  setup  boundary  conditions, 
prepare  input  required  by  the  solver,  perform  domain  decomposition  for  parallel  processing,  and 
execute  the  solver.  Our  design  of  a  common  platform  for  multidisciplinary  analysis  is  based  on 
our  experience  in  developing  TEMPUS. 
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Figure  17.  TEMPUS  environment  for  wide  band  CEM  applications 


At  the  core  of  the  TEMPUS  environment,  there  is  a  collection  of  technologies  that  represent  the 
strength  of  TEMPUS.  Some  of  them  are  listed  here: 

1.  A  higher  order  discontinuous  Galerkin  method  that  is  currently  operational  up  to  10* 
order  for  spatial  representation  of  electric  and  magnetic  fields  inside  each  computational 
finite-volume  cell. 

2.  A  hybrid  unstructured  grid  framework  for  modeling  complete  targets  including  engine 
inlets  with  stages  of  fan  blades 

3.  Modeling  of  general  materials  such  as  lossy/lossless  dielectric,  resistive  cards,  impedance 
layers  and  dispersive  media 

4.  Highly  scalable  parallel  code  architecture  that  is  transportable  across  different  platforms. 


11.3  Higher  Order  Schemes  for  Multidisciplinary  Analysis 

Our  choice  for  higher-order  multidisciplinary  analysis  is  the  DG  scheme.  This  choice  is  based 
not  only  on  our  extraordinary  success  in  developing  a  DG  based  higher-order  solver  for  CEM 
that  clearly  demonstrates  the  ability  of  a  higher-order  DG  scheme  to  produce  very  accurate 
solutions  with  appreciable  decrease  in  turnaround  time,  but  also  on  the  fact  that  DG  is  the  most 
appropriate  choice  for  multidisciplinary  analysis  based  on  the  discussion  presented  in  the  next 
three  paragraphs. 
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One  of  the  earliest  attempts  at  developing  highly  accurate  numerical  schemes  employed  a 
spectral  method.  In  a  spectral  method  the  dependent  variables  Q  are  expanded  in  terms  of  a 

basis  function  {v^{x,y,z),i  =  l,nj  as  Q=  .  The  solution  process  solves  for  the  coefficients 

i=l 

ffi  =  \,n.  This  scheme  requires  the  solution  to  be  continuous  across  the  computational  domain 
and  performs  exceedingly  well  for  problems  that  satisfy  the  requirement.  Unfortunately  a  large 
category  of  problems  that  we  encounter  in  real  life  does  not  satisfy  such  a  requirement.  A 
natural  extension  to  the  spectral  method  that  is  less  restrictive  is  to  employ  a  local  basis  function 
rather  than  a  global  one.  The  DG  scheme  does  just  that.  Use  of  a  local  basis  function  permits 
the  solution  to  be  discontinuous  at  element  boundaries  and  thus  enables  application  of  the 
scheme  to  problems  wherein  accurate  solutions  may  be  obtained  by  resolving  discontinuities  at 
element  interfaces.  Within  each  element  the  solution  is  assumed  to  be  continuous.  This 
restriction  is  obviously  much  less  severe  than  the  restriction  imposed  by  a  spectral  method 
wherein  the  solution  is  required  to  be  continuous  across  the  whole  computational  domain. 

Essentially  non-oscillatory  (ENO)  schemes  and  the  weighted  ENO  (WENO)  schemes  also 
employ  local  basis  functions.  The  main  differences  between  the  ENO  schemes  and  the  DG 
schemes  are: 

1 .  The  ENO  schemes  solve  for  the  dependent  variables  while  the  DG  schemes  solve  for  the 
coefficients  QU  =  l,n 

2.  In  an  ENO  scheme  the  coefficients  Q'^i  =  \,n  are  computed  by  an  appropriate  fit  of  the 
data  in  the  neighborhood  of  an  element  while  in  a  DG  scheme  a  weak  form  of  the 
solution  is  employed  in  the  computation  of  the  coefficients. 

Item  1  above  implies  that  the  DG  scheme  is  solving  for  more  unknowns  than  the  ENO  schemes. 
Item  2  implies  that  the  ENO  schemes  employ  a  much  larger  numerical  stencil  than  the  DG 
schemes.  As  advantageous  as  it  is  to  solve  for  less  number  of  unknowns,  the  use  of  a  much 
larger  stencil  in  ENO  schemes  renders  the  schemes  almost  unusable  due  to  severe  stability 
limitations. 

One  of  the  attractive  features  of  the  DG  schemes  is  the  compactness  of  the  stencil  employed. 
That  is,  the  equations  that  determine  the  coefficients  =  \,n  depend  only  on  the  data  in  the 
immediate  neighborhood  of  a  stencil.  The  FADE  type  of  finite  difference  schemes  also  has  this 
property  of  compactness.  Unfortunately  such  schemes  do  not  extend  easily  to  unstructured 
grids.  The  DG  schemes  extend  quite  naturally  to  discretizing  the  computational  domain  using 
elements  of  arbitrary  shapes  and  as  such  have  a  clear  advantage  over  the  FADE  schemes.  The 
advantages  that  unstructured  grids  have  over  structured  grids  are  so  powerful  that  it  is  almost 
imperative  that  we  choose  a  scheme  that  extends  easily  to  unstructured  grids. 

Based  on  the  above  discussion,  the  higher-order  multidisciplinary  platform  that  we  will  develop 
will  employ  a  DG  scheme  and  unstructured  grids. 
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11.4  Computational  Environment  for  Multidisciplinary  Analysis 

There  are  two  prineipal  methods  to  integrate  multiple  physical  phenomena  into  a  single 
computational  environment.  They  are: 

(a)  Tight  (monolithic)  coupling,  in  which  the  various  physical  modules  are  integrated  into 
one  large  matrix  system  using  a  single  numerical  strategy  and  inverted  in  bulk,  and, 

(b)  Loose  (iterative)  coupling,  in  which  the  various  physical  models  are  solved  individually 
and  global  iterations  are  performed  in  order  to  converge  them  in  a  given  time  step 

Most  of  the  multidisciplinary  analyses  currently  carried  out  employ  loosely  coupled  algorithms. 
For  example,  in  the  analysis  of  the  interaction  between  aerodynamic  and  elastic  forces, 
aerodynamic  loads  are  computed  first  and  the  resulting  elastic  response  is  then  evaluated  to 
obtain  a  modified  shape  which  is  used  to  recompute  the  aerodynamic  loads.  This  process  is 
continued  until  a  convergence  criterion  is  satisfied.  As  such  a  loosely  coupled  algorithm  only 
requires  an  interface  between  different  solvers  employed  in  the  simulation  of  different  physical 
phenomena.  This  may  be  quite  easily  achieved  by  developing  a  modular  environment  and  the 
required  interfaces.  Our  goal  is  to  go  beyond  the  development  of  such  an  environment  and 
develop  strongly  coupled  algorithms  wherever  possible.  Magnetogasdynamics,  multiphase,  and 
free  surface  flows  are  examples  wherein  a  strongly  coupled  algorithm  may  be  developed.  Fig.  18 
shows  a  sampling  of  multidisciplinary  problems  that  may  be  handled  by  the  proposed 
environment. 

Although  a  loosely  coupled  multidisciplinary  analysis  may  appear  attractive  for  its  simplicity  in 
implementation,  freedom  in  discretization  procedures,  and  ease  in  extending  to  different 
phenomena,  a  tightly  coupled  analysis  is  certainly  desirable  since  it  is  likely  to  be  more  accurate, 
more  stable,  and  will  permit  use  of  larger  time-steps  for  time-accurate  simulations.  Our  choice  is 
a  combination  of  two.  We  propose  to  employ  a  tightly  coupled  strategy  wherever  possible.  For 
example,  a  multi-species  flow  with  chemical  reaction  may  be  solved  in  a  coupled  manner  while 
the  analysis  of  the  interaction  between  aerodynamic  and  structural  forces  will  require  a  loosely 
coupled  strategy. 

One  of  the  challenges  that  we  need  to  address  is  related  to  optimizing  the  resource  requirements. 
In  a  higher-order  DG  scheme,  the  number  of  dependent  variables  increases  with  the  order  of  the 
scheme.  Since  the  computing  resources  required  are  usually  nonlinear  functions  of  the  number 
of  dependent  variables,  it  is  quite  easy  to  stress  the  available  resources  beyond  their  limits.  The 
main  task  is  to  ensure  that  we  employ  an  adaptive  order  scheme  so  that  we  use  only  as  high  an 
order  of  scheme  for  an  element  as  is  required  by  the  solution  and  not  any  higher. 

The  second  challenge  is  developing  an  efficient  strategy  for  parallel  execution  of 
multidisciplinary  simulation.  In  the  case  of  a  tightly  coupled  strategy,  this  is  not  a  major  issue 
since  the  same  grid  is  employed  in  solving  different  governing  equations.  Any  standard  domain 
decomposition  module  such  as  METIS  may  be  employed  to  ensure  efficient  parallel 
performance.  In  the  case  of  a  loosely  coupled  system,  we  first  need  to  evaluate  the  CPU  time 
required  by  different  solvers  and  based  on  the  grid  employed  and  the  time  required  by  the  solver 
we  need  to  develop  an  optimal  distribution  of  work  among  available  processors.  In  order  to 
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ensure  maximum  efficiency,  we  may  have  to  perform  different  number  of  iterations  for  different 
solvers  before  data  is  exchanged  between  them. 
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Figure  18.  Some  Multiphysics  Problems  of  Interest 


Figs.  19  and  20  show  a  sample  strategy  to  parallelize  an  incompressible  free  surface  MHD  flow 
simulation.  Various  physical  elements  are  individually  converged  to  an  acceptable  tolerance. 
Different  physical  phenomena  may  require  different  grid  resolution,  thus  requiring  a  strategy  to 
exchange  information  across  grid  levels.  As  discussed  in  section  11.1,  use  of  DG  as  a  higher- 
order  scheme  simplifies  exchange  of  information  between  different  grids  and  the  only  module 
that  needs  to  be  developed  is  an  octree  based  technique  for  locating  the  element  in  one  grid  that 
contain  a  given  point  from  another  grid. 
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Figure  19.  Pre-processing  setup  for  HIMAG:  HyPerComp  Incompressible 
MHD  solver  for  Arbitrary  Geometries 


time  level  n 


Momentum  Solver  (getdu) 


time  level  n+1 


Figure  20.  Parallelization  strategy  for  the  incompressible  MHD  solver 
HIMAG  -  the  8  s  represent  tolerances  of  individual  solvers 
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12.0  Suggestions  for  future  work 

Our  proposal  for  the  continuation  of  the  work  done  under  Phase  I  involves  the  following  3 
categories: 

1.  Construction  of  a  general  purpose  three-dimensional  code  environment  for  multiphysics 
simulation.  The  starting  point  for  this  task  is  the  existing  TEMPUS  CEM  environment. 

2.  Implementation  and  enhancement  of  high  order  algorithms  addressing  multiphysics 

3.  Demonstration  of  the  high  order  platform  for  selected  candidate  problems  requiring 
multiphysics  simulations 

Detailed  list  of  subtasks  to  be  done  under  each  of  the  tasks  listed  above  is  given  below. 

T-1:  Multidisciplinary  Environment  Development  Tasks 

T-1.1:  Starting  with  the  TEMPUS  CEM  environment,  develop  and  incorporate  data  structure  for 
modeling  multiphysics  simulations  and  optimize  memory  resource  requirements 
T-1. 2:  Develop  GUIs  for  problem  set  up,  boundary  conditions  and  run  set  up  for  various  physics 
T-I.3:  Develop  interfaces  including  interpolation  routines  for  data  communications  between 
different  solvers  for  different  pysics. 

T-1. 4:  Develop  appropriate  parallel  MPI  routines  for  multiphysics  data  communication  across 
nodal  boundaries 

1 .4. 1  Parallel  data  structure  for  moving  mesh 

1 .4.2  Load  balancing  for  moving/adaptive  mesh 

T-1. 5:  Add  to  the  GUI  a  database  for  relevant  physical  properties 

T-I.6:  Develop  makefile/compiler  options  for  portability  across  different  parallel  hardware 

T-2;  High  Order  Algorithm  Enhancement  Tasks 

T-2.1:  Develop  orthonormal  basis  function  set  for  at  least  up  to  10*  order  for  multiphysics 
T-2.2:  Special  basis  functions  for  singular  regions  (e.g.  sharp  leading  edge  regions) 

T-2. 3:  High  order  boundary  condition  procedures  (e.g.  nonreflecting  outer  boundary  conditions) 
T-2.4:  Develop  linearization  procedures  for  advection  terms 
T-2. 5:  Develop  Riemann  flux  routines  for  different  physics 
T-2. 6:  Multi-order  and  implicit  time  step  procedures 

T-2. 7:  Limiters/filters/pre-conditioners/penalty  functions  and  other  stabilizing/code  acceleration 
techniques 

T-2. 8:  Incorporation  of  high  order  geometry  representation  and  nonorthogonal  grid  effects 
T-2. 9:  Coupling  procedures  for  elliptic,  parabolic  and  hyperbolic  regions/aspects  of 
multiphysics 

T-2.10:High  order  level  set  procedures  for  free  surface  tracking 

T-3;  Tasks  on  Multiphysics  Simulation/Demonstration 

This  set  of  tasks  will  complement  ongoing  research  efforts  and  will  synchronize  with  the  time 
lines  in  which  model  development  from  those  projects  will  supplement  and  validate  present 
results.  Unit  problems  in  MHD,  turbulent  fluid  flow  ,  acoustics  and  aerothermal  flows  will  be 
considered.  The  first  two  problems  below  will  be  used  as  demonstration  cases  on  large  scale 
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parallel  LINUX  cluster.  Specifics  of  the  later  tasks  will  be  determined  jointly  with  the  program 
monitor. 

T-3. 1 :  Euler/Navier-Stokes/Maxwell  equations  for  MGD  with  chemical  kinetics 
Hypersonic  inlet  flow  with  viscosity,  (RANS)  turbulence  and  MHD  terms. 

T-3. 2:  Incompressible  flow  with  MHD  and  a  free  surface 

Free  surface  flow  of  liquid  metal  past  a  cylinder  in  the  presence  of  a  magnetic  field 
T-3. 3:  Demonstration  of  high  order  DG  for  a  candidate  turbulent  flow  problem 
T-3. 4:  Demonstration  of  a  flexible  structure  with  active  feedback  control 
T-3. 5:  Study  of  high  order  absorbing  outer  boundary  conditions  for  an  acoustics  problem 
T-3. 6:  Application  of  high  order  DG  for  accurate  simulation  of  aerothermal  prediction 

T-4;  Analysis 

In  addition  to  all  of  the  above,  we  intend  to  perform  preliminary  analysis  of  the  DG  scheme  as 
we  develop  it,  and  study  its  stability  and  convergence  properties. 
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Appendix  1 

Flowchart  starting  from  main.c 
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Read  in  input 

Read  in  grid 
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Read  in  partition  boundary 
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